************************************************************************
* Project: The Press-Safety Paradox of Democracies in Media Systems,   *
*          Regime-type Durability and Journalist Killings              *
* Author: Jonathan A. Solis                                            *
* Task: Replication Materials for Foreign Policy Analysis article      *
* Stata version: IC 15.1                                               *
* Date: 7/24/2019                                                      *
************************************************************************
*Creat log file
*log using Solis_fpa_Stata

*Load Data
clear 
cd "[working directory here]"
import delimited "Solis_fpa_july2019"
set more off

*Prepare dataset for analysis*

*sort 
sort country year

*n. log regime durability and other variables
gen seq_ln=ln(seq1)
gen gdppc_ln=ln(gdppc)
gen gdp_ln=ln(gdp)

*create dummy years (will need separate variables for rare events logit/figures)
forvalues i=1992/2016{

gen y`i'=1 if year==`i'
replace y`i'=0 if y`i'==.
}
*

*Create binary journalist killed variable
gen logit=1 if confirmed!=0 & confirmed!=.
replace logit=0 if logit==.

*Create journalist killed ordinal variables
  
*Ordinal 1 (not killed, one killed, 2 or more killed)
  gen ord1=.
  replace ord1=0 if conf ==0 
  replace ord1=1 if conf==1
  replace ord1=2 if conf>=2

 *Ordinal 2 (not killed, 1-9 killed, 10 or more killed)
  gen ord2=.
  replace ord2=0 if conf ==0
  replace ord2=1 if inrange(conf,1,9)
  replace ord2=2 if conf>=10

*-----------------------*
************************* 
***MANUSCRIPT ANALYSIS***
************************* 
*-----------------------*

*********  
*Table 1*
*********
  
*With original Asal et al. (2018) indicators; 1992-2011*
  
*Model 1: Logit
logit logit seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*Store estimates
est sto m1

*Percentage of zeros? (how many country-years without a journalist killed)
predict resid
replace resid=. if year>=2012
sum resid
gen missing1=.
replace missing1 = 1 if resid!=.
replace missing1 = 0 if resid==.
tab missing1
*about 89% 0's

*AIC
estat ic
*AIC=1224.296

*Model 2: Rare Events Logit
relogit logit seq_ln polity qog physint speech intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 ///
  y2002 y2003 y2004 y2005 y2006 y2007 y2008 y2009 y2010, cluster(ccode) 

*Store estimates
est sto m2

*Unable to get AIC from RE Logit model

*Model 3: Ordinal 1
ologit ord1 seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*Store estimates
est sto m3

*AIC
estat ic
*AIC=1562.543

*Model 4: Ordinal 2
ologit ord2 seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*Store estimates
est sto m4

*AIC
estat ic
*AIC=1258.791
  
*Model 5: NBREG
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*Store estimates
est sto m5

*AIC
estat ic
*AIC=1945.725

*Model 6: ZINB
zinb confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, inflate(pop_ln) cluster(ccode) nolog

*Store estimates
est sto m6

*AIC
estat ic
*AIC=1949.184

****Create TABLE 1 for LaTex (basic table; I make some changes by hand once generated)****
esttab m1 m2 m3 m4 m5 m6 using table1.tex, replace se aic obslast r2 ///
mtitle("Logit" "RE Logit" "Ordinal 1" "Ordinal 2" "NBREG" "ZINB" ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" qog ///
"Quality of Govt." physint "Physical Integrity" speech "Freedom of Speech" ///
 intensity2 "Armed Conflict" info "Information Flows" pop_ln "Population (ln)") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Generate Fig. 2

*run model 1 w/ year variables
quietly nbreg confir seq_ln polity qog physint speech intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 ///
  y2002 y2003 y2004 y2005 y2006 y2007 y2008 y2009 y2010, cluster(ccode)

*create figure
prgen seq_ln, from (0) to (5.5) generate(durat) rest(mean) ci 
label variable duratp1 "Probability of Journalist Killing" 
label variable duratx "Regime Duration"  
label variable duratp1lb "95% lower limit" 
label variable duratp1ub "95% upper limit"

twoway (connected duratp1 duratx,lpattern(solid) lwidth(thin) msymbol(none) ///
yaxis(1)) (connected duratp1lb duratx, lpattern(dash) lwidth(thin) msymbol(none) ///
yaxis(1)) (connected  duratp1ub duratx, lpattern(dash) lwidth(thin) msymbol(none) yaxis(1) ///
ylab(,nogrid))

*Note: I further stylize figure for manuscript
  
*********
*Table 2*
*********

*Expanded models w/ new variables (1992-2014)

*NBREG

*Model 7: All regime types (worldwide)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog 

*Identify observation in full sample (will need for appendix analysis) 
gen byte full=e(sample)

*Store estimates
est sto t2m1

*AIC
estat ic
*AIC=2860.384 

*IRR=.7329813 (seq_ln)

************************************
*Incidence Rate Ratio (IRR) Results*
************************************

*Add the 'irr' after comma to reproduce Incidence Rate Ratio (IRR) results. For example,
* for model 7:
*
*nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
*  y199* y2*, irr cluster(ccode) nolog 
*IRR=.7329813 (seq_ln)

*Model 8: Only autocracies
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0, cluster(ccode) nolog

*Identify observation in autocracies sample (will need for appendix analysis) 
gen byte auto=e(sample)

*Store estimates
est sto t2m2

*AIC
estat ic
*AIC=419.0924

*IRR=.6551315 (seq_ln)

*Model 9: Only anocracies
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1, cluster(ccode) nolog

*Identify observation in anocracies sample (will need for appendix analysis) 
gen byte ano=e(sample)

*Store estimates
est sto t2m3

*AIC
estat ic
*AIC=1039.686

*IRR=.5550174 (seq_ln)

*Model 10: Only democracies
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2, cluster(ccode) nolog

*Identify observation in anocracies sample (will need for appendix analysis) 
gen byte demo=e(sample)

*Store estimates
est sto t2m4

*AIC
estat ic
*AIC=1346.538

*IRR=1.072929 (seq_ln)
*IRR=.7273303 (polity level)

****Create TABLE 2 for LaTex (basic table; I make some changes by hand once generated)****
esttab t2m1 t2m2 t2m3 t2m4 using table2.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Generate Fig. 3

*run model 7
quietly nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 y2002 y2003 y2004 ///
  y2005 y2006 y2007 y2008 y2009 y2010 y2011 y2012 y2013, cluster(ccode) nolog 

*drop durt* variables if you have not cleared from fig. 2 
drop durat*

*create figure
prgen seq_ln, from (0) to (5.5) generate(durat) rest(mean) ci 
label variable duratp1 "Probability of Journalist Killing" 
label variable duratx "Regime-type Duration"  
label variable duratp1lb "95% lower limit" 
label variable duratp1ub "95% upper limit"

twoway (connected duratp1 duratx, lpattern(solid) lwidth(thin) msymbol(none) ///
yaxis(1)) (connected duratp1lb duratx, lpattern(dash) lwidth(thin) msymbol(none) ///
yaxis(1)) (connected  duratp1ub duratx, lpattern(dash) lwidth(thin) msymbol(none) yaxis(1) ///
ylab(,nogrid))

*Note: I further stylize figure for manuscript

********************************************************************
* SEE R SCRIPT FOR DURATION ANALYSIS LABELED `Solis_fpa_aug2019.R' *
********************************************************************

************************************
***********APPENDIX*****************
************************************

***********************************
***********************************
** SECTION A: Summary Statistics **
***********************************
***********************************

*Create variable for unconfirmed journalist killed +1: Zero Inflated NBREG does
* not run otherwise
gen fix=unconfirmed+1
corr fix unconfirmed

*             |      fix unconf~d
*-------------+------------------
*         fix |   1.0000
* unconfirmed |   1.0000   1.0000

*Note, variables full, auto, ano, and demo created above; indicate their
* respective samples

*Main Analysis (A.1)
sum confirmed seq_ln polity2 qog physint speech intensity2 info pop_ln ///
 public_cor express_vd physical_vd 
 
*Robustness checks (A.1)
sum fix gmf fhfp gdp_ln gdppc_ln gdppc_cng hom_count hom_rate 
 
*Full (A.2)
sum confirmed seq_ln seq1 polity2  intensity2 info pop_ln ///
 public_cor express_vd physical_vd if full==1
 
*Autocracies (A.2)
sum confirmed seq_ln seq1 polity2 intensity2 info pop_ln ///
 public_cor express_vd physical_vd if auto==1

*Anocracies (A.2)
sum confirmed seq_ln seq1 polity2  intensity2 info pop_ln ///
 public_cor express_vd physical_vd if ano==1

*Democracies (A.2)
sum confirmed seq_ln seq1 polity2 intensity2 info pop_ln ///
 public_cor express_vd physical_vd if demo==1

*Correlation Matrix (A.3)
pwcorr qog  public_cor 
pwcorr physint physical_vd 
pwcorr speech express_vd  

*************************************
*************************************
** SECTION B: DESCRIPTIVE ANALYSIS **
*************************************
*************************************

*Data and code for descriptive analysis available upon request.

*Figure out countries with no polity score in sample (B.2 intro)
****Countries that have missing polity
list country year if polity==.  & confirmed>=1 & inrange(year,1992,2014)
	
*      +--------------------+
*      |     country   year |
*      |--------------------|
*  10. | Afghanistan   2001 |
*  15. | Afghanistan   2006 |
*  16. | Afghanistan   2007 |
*  17. | Afghanistan   2008 |
*  18. | Afghanistan   2009 |
*      |--------------------|
*  19. | Afghanistan   2010 |
*  20. | Afghanistan   2011 |
*1657. |        Iraq   2003 |
*1658. |        Iraq   2004 |
*1659. |        Iraq   2005 |
*      |--------------------|
*1660. |        Iraq   2006 |
*1661. |        Iraq   2007 |
*1662. |        Iraq   2008 |
*1663. |        Iraq   2009 |
*1998. |     Lebanon   1992 |
*      |--------------------|
*1999. |     Lebanon   1993 |
*2005. |     Lebanon   1999 |
*3198. |     Somalia   2011 |
*      +--------------------+
			

*Number we lose from these
tab confirmed if polity==.  & confirmed>=1. & inrange(year,1992,2014)

*  confirmed  Freq.	  confirmed*Freq.
      *1	    2	        2
      *2	    7	        14
      *3	    1	        3
      *4	    1	        4
      *9	    1	        9
      *11	    1	        11
      *14	    1 	        14
      *23	    1 	        23
      *24	    1 	        24
      *32	    2 	        64
	                  *------------*
	                    *total=168
						
***************************************************
***************************************************
** SECTION C: MODEL'S CONTROL VARIABLES DECISOIN **
***************************************************
***************************************************

*Generate Table 5*

*Model 1: Global (NBREG)
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*Store estimates and create variable to identify sample
est sto rra_full
gen byte full_asalb=e(sample)

*AIC
estat ic
*AIC=1945.725

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2* if durable2==0, cluster(ccode) nolog

*Store estimates and create variable to identify sample
est sto rra_auto
gen byte auto_asalb=e(sample)

*AIC
estat ic
*AIC=253.645

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2* if durable2==1, cluster(ccode) nolog

*Store estimates and create variable to identify sample
est sto rra_ano
gen byte ano_asalb=e(sample)

*AIC
estat ic
*AIC=587.702

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2* if durable2==2, cluster(ccode) nolog

*Store estimates and create variable to identify sample
est sto rra_demo
gen byte demo_asalb=e(sample)

*AIC
estat ic
*AIC=1120.101

****Create TABLE 5 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab rra_full rra_auto rra_ano rra_demo using apx_revres1.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" qog ///
"Quality of Govt." physint "Physical Integrity" speech "Freedom of Speech" ///
 intensity2 "Armed Conflict" info "Information Flows" pop_ln "Population (ln)") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Generate Table 6*

*Model 1: NBREG
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if inrange(year,1992,2011), cluster(ccode) nolog 

*Store estimates and create variable to identify sample  
est sto rr_full
gen byte full_solisb=e(sample)

*AIC
estat ic
*AIC=2425.738 

*Model 2: NBREG
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0 & inrange(year,1992,2011), cluster(ccode) nolog 

*Store estimates and create variable to identify sample
est sto rr_auto
gen byte auto_solisb=e(sample)

*AIC
estat ic
*AIC=365.3766

*Model 3: NBREG
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1 & inrange(year,1992,2011), cluster(ccode) nolog 

*Store estimates and create variable to identify sample
est sto rr_ano
gen byte ano_solisb=e(sample)

*AIC
estat ic
*AIC=853.746

*Model 4: NBREG
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2 & inrange(year,1992,2011), cluster(ccode) nolog

*Store estimates and create variable to identify sample
est sto rr_demo
gen byte demo_solisb=e(sample)

*AIC
estat ic
*AIC=1163.697

****Create TABLE 6 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab rr_full rr_auto rr_ano rr_demo using rr2011_solis.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" gmf "Global Media Freedom" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" polity2 "Polity" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 gdp_ln "GDP (ln)" gdppc_ln "GDP p/c (ln)" gdppc_cng "$\Delta$ GDP p/c" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Generate data for Table 7*

*Find data for percentage less of observations
*Instructions: take samples sizes from Appendix tables 5 and 6, then compare

*Global samples:
* asal controls = 2,491 
* vdem controls = 3,114
display 3114 - 2491
* = 623
display (623/3114)*100
*~20%

*Auto samples:
* asal controls = 391 
* vdem controls = 546
display 546 - 391
* = 155
display (155/546)*100
*~28%

*Ano samples:
* asal controls = 643 
* vdem controls = 920
display 920 - 643
* = 277
display (277/920)*100
*~30%

*Demo samples:
* asal controls = 1457 
* vdem controls = 1648
display 1648 - 1457
* = 191
display (191/1648)*100
*~11%

*Find data for percentage less of journalists killed in missing observations

*Rerun models for appendix table 5 and 6

*appendix table 5
nbreg confirmed seq_ln polity qog physint speech intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog
  
*Create variable to identify sample
gen byte full_asal=e(sample)

*appendix table 6
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if inrange(year,1992,2011), cluster(ccode) nolog 

*Create variable to identify sample
gen byte full_solis=e(sample)

*Instructions: tabulate journalists killed, multiply `confirmed' columns and `Freq.'
* column. Then add these results.

*Number of journalists killed in table 5 samples 
tab confirmed if full_asal==1 & inrange(year,1992,2011)
*587
tab confirmed if full_asal==1 & inrange(year,1992,2011) & durable2==0 
*58
tab confirmed if full_asal==1 & inrange(year,1992,2011) & durable2==1 
*181
tab confirmed if full_asal==1 & inrange(year,1992,2011) & durable2==2 
*348

*Number of journalists killed in table 6 samples 
tab confirmed if full_solis==1 & inrange(year,1992,2011)
*731
tab confirmed if full_solis==1 & inrange(year,1992,2011) & durable2==0 
*101
tab confirmed if full_solis==1 & inrange(year,1992,2011) & durable2==1 
*271
tab confirmed if full_solis==1 & inrange(year,1992,2011) & durable2==2 
*359

*Generate Table 8*

*Generate variable for missingness
gen missing=.
replace missing =1 if full_solis==1 & full_asal!=1
replace missing=0 if missing==.

*Model likelihood of missingness w/polity
logit missing polity   ///
  y199* y2* if full_solis==1 & inrange(year,1992,2011),  cluster(ccode) nolog 

*Store estimates
est sto missing1

*Predict probability of missingness
predict pr_a if e(sample)
*predict idxhat if e(sample), xb
*generate phat2 = exp(idxhat)/(1+exp(idxhat))

*Model likelihood of missingness w/polity and other controls
logit missing polity pop_ln gdppc_ln intensity2  ///
  y199* y2* if full_solis==1 & inrange(year,1992,2011), cluster(ccode) nolog 
est sto missing2

*Predict probability of missingness
predict pr_b if e(sample)
*predict idxhat if e(sample), xb
*generate phat2 = exp(idxhat)/(1+exp(idxhat))

****Create TABLE 8 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab missing1 missing2 using missingness.tex, replace se aic obslast r2 ///
coeflabel( polity2 "Polity" pop_ln "Population (ln)" gdppc_ln "GDP p/c (ln)" ///
 intensity2 "Conflict") ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Generate data Table 9*

*Get predicted probabilities for models in table 8 by regime type

*Model 1, table 8
*Avg by regime type
* -> Get the mean from these commands
sum pr_a if durable2==0
sum pr_a if durable2==1
sum pr_a if durable2==2

*Model 2, table 8
*Avg by regime type
* -> Get the mean from these commands
sum pr_b if durable2==0
sum pr_b if durable2==1
sum pr_b if durable2==2

*******************************************************************
*******************************************************************
** SECTION D: COX REGRESSION MODEL MULTIPLE EVENT SPECIFICATIONS **
*******************************************************************
*******************************************************************

* SEE R SCRIPT FOR SURVIVAL ANALYSIS LABELED `Solis_fpa_may2019.R' 

***********************************************
***********************************************
** SECTION E: ALTERNATE MODEL SPECIFICATIONS **
***********************************************
***********************************************

*Table 11: Logit  

*Model 1: Global
logit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* , cluster(ccode) nolog

*store estimates
est sto log1

*AIC  
estat ic
*AIC=2860.384

*Model 2: Auto
logit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates
est sto log2

*AIC  
estat ic
*AIC=419.092

*Model 3: Ano
logit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1, cluster(ccode) nolog

*store estimates
est sto log3

*AIC
estat ic
*AIC=1039.686

*Model 4: Demo
logit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2, cluster(ccode) nolog

*store estimates
est sto log4

*AIC
estat ic
*AIC=1346.538

****Create TABLE 11 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab log1 log2 log3 log4 using apx_logit.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" unconfirmed "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 12: Rare Events Logit

*Model 1: Global
relogit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 y2002 y2003 y2004 ///
  y2005 y2006 y2007 y2008 y2009 y2010 y2011 y2012 y2013, cluster(ccode)

*store estimates
est sto rel1

*Model 2: Auto
relogit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 y2002 y2003 y2004 ///
  y2005 y2006 y2007  y2009 y2010 y2011 y2012 y2013 if durable2==0, cluster(ccode)

*store estimates
est sto rel2

*Model 3: Ano
relogit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 y2002 y2003 y2004 ///
  y2005 y2006 y2007 y2008 y2009 y2010 y2011 y2012 y2013 if durable2==1, cluster(ccode)

*store estimates  
est sto rel3

*Model 4: Demo
relogit logit seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y1992 y1993 y1994 y1995 y1996 y1997 y1998 y1999 y2000 y2001 y2002 y2003 y2004 ///
  y2005 y2006 y2007 y2008 y2009 y2010 y2011 y2012 y2013 if durable2==2, cluster(ccode)

*store estimates
est sto rel4

****Create TABLE 12 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab rel1 rel2 rel3 rel4 using apx_rel.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" unconfirmed "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz
  
*Table 13: Ordinal 1

*Model 1: Global
ologit ord1 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog
 
*store estimates
est sto ordv1_1

*AIC
estat ic
*AIC=2198.96

*Model 2: Auto
ologit ord1 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates
est sto ordv1_2

*AIC
estat ic
*AIC=339.4809

*Model 3: Ano
ologit ord1 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1, cluster(ccode) nolog

*store estimates
est sto ordv1_3

*AIC
estat ic
*AIC=821.8949

*Model 4: Demo
ologit ord1 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2, cluster(ccode) nolog

*store estimates
est sto ordv1_4

*AIC
estat ic
*AIC=1047.007

****Create TABLE 13 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab ordv1_1 ordv1_2 ordv1_3 ordv1_4 using apx_ord1.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" unconfirmed "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 14: Ordinal 2

*Model 1: Global
ologit ord2 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2*, cluster(ccode) nolog

*store estimates
est sto ordv2_1

*AIC
estat ic
*AIC=1820.844

*Model 2: Auto
ologit ord2 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates
est sto ordv2_2

*AIC
estat ic
*AIC=313.3774

*Model 3: Ano
ologit ord2 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1, cluster(ccode) nolog

*store estimates
est sto ordv2_3

*AIC
estat ic
*AIC=691.2462

*Model 4: Demo
ologit ord2 seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2, cluster(ccode) nolog

*store estimates 
est sto ordv2_4

*AIC
estat ic
*AIC=810.4865

****Create TABLE 14 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab ordv2_1 ordv2_2 ordv2_3 ordv2_4 using apx_ord2.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" unconfirmed "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 15: ZINB

*Model 1: Global
zinb confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* , inflate(pop_ln fix) cluster(ccode) nolog

*store estimates 
est sto zinb1

*AIC
estat ic
*AIC=2810.225

*Model 2: Auto
zinb confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==0, inflate(pop_ln fix) cluster(ccode)

*store estimates 
est sto zinb2

*AIC
estat ic
*AIC=394.562

*Model 3: Ano
zinb confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==1, inflate(pop_ln fix) cluster(ccode) nolog

*store estimates 
est sto zinb3

*AIC
estat ic
*AIC=1033.192

*Model 4: Demo
zinb confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln ///
  y199* y2* if durable2==2, inflate(pop_ln fix) cluster(ccode) nolog

*store estimates 
est sto zinb4

*AIC
estat ic
*AIC=1334.986

****Create TABLE 15 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab zinb1 zinb2 zinb3 zinb4 using apx_zinb.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 16: LDV Logit

*Set for time series to activate lag command (l.) 
tsset ccode year

*Model 1: Full
nbreg confirmed l.confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln   ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates 
est sto lag1

*AIC
estat ic
*AIC=2694.867

*Model 2: Auto
*nbreg confirmed l.confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  ///
  *  y199* y2* if durable2==0, cluster(ccode) 

  *maximum likelihood estimator is not concave and fails to produce a
* maximum likelihood estimate	

*Model 3: Ano
nbreg confirmed l.confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates 
est sto lag3

*AIC
estat ic
*AIC=969.388

*Model 4: Demo
nbreg confirmed l.confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates 
est sto lag4

*AIC
estat ic
*AIC=1304.015

****Create TABLE 16 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab lag1  lag3 lag4 using apx_ldv.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(L.confirmed "Lagged DV" seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed") ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

**********************************
**********************************
** SECTION F: ROBUSTNESS CHECKS **
**********************************
**********************************

*Table 17: w/ Global Media Freedom (GMF) covariate

*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  gmf ///
  y199* y2*, cluster(ccode)  nolog 

*store estimates 
est sto gmf1

*AIC
estat ic
*AIC=2855.827

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gmf ///
    y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates 	
est sto gmf2

*AIC
estat ic
*AIC=421.008

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  gmf  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates 	
est sto gmf3

*AIC
estat ic
*AIC=1037.372

*Model 4: Demo (NBREG)
  nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gmf  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates 	
est sto gmf4

*AIC
estat ic
*AIC=1347.419

****Create TABLE 17 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab gmf1 gmf2 gmf3 gmf4 using apx_gmf.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" gmf "Global Media Freedom" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" polity2 "Polity" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 gdp_ln "GDP (ln)" gdppc_ln "GDP p/c (ln)" gdppc_cng "$\Delta$ GDP p/c" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 18: w/ Freedom House's Press Freedom in the World covariate

*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln fhfp ///
  y199* y2*, cluster(ccode)  nolog 

*store estimates 	
est sto fh1

*AIC
estat ic
*AIC=2684.999

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln fhfp  ///
   y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates 	
est sto fh2

*AIC
estat ic
*AIC=381.058

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  fhfp   ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates 	
est sto fh3

*AIC
estat ic
*AIC=965.771

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln fhfp   ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates 	
est sto fh4

*AIC
estat ic
*AIC=1307.2

****Create TABLE 18 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab fh1 fh2 fh3 fh4 using apx_fh.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" fhfp "Press Freedom, Freedom House" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" polity2 "Polity" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 gdp_ln "GDP (ln)" gdppc_ln "GDP p/c (ln)" gdppc_cng "$\Delta$ GDP p/c" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 19: w/ GDP covariate

*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  gdp_ln ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates 	
est sto gdp1

*AIC
estat ic
*AIC=2643.251

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info  pop_ln gdp_ln ///
    y199* y2*   if durable2==0, cluster(ccode)

*store estimates 	
est sto gdp2

*AIC
estat ic
*AIC=339.576 (but changes, for some reason)

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info  pop_ln gdp_ln  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto gdp3

*AIC
estat ic
*AIC=910.8602

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gdp_ln ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates
est sto gdp4

*AIC
estat ic
*AIC=1329.935

****Create TABLE 19 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab gdp1 gdp2 gdp3 gdp4 using apx_gdp.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 gdp_ln "GDP (ln)" gdppc_ln "GDP p/c (ln)" gdppc_cng "$\Delta$ GDP p/c" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 20: w/ GDP per capita covariate

*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gdppc_ln gdppc_cng ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates
est sto gdppc1

*AIC
estat ic
*AIC=2571.473

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info  gdppc_ln gdppc_cng ///
    y199* y2*   if durable2==0, cluster(ccode)

*store estimates
est sto gdppc2

*AIC
estat ic
*AIC=330.179 (but changes, for some reason)

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gdppc_ln gdppc_cng   ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto gdppc3

*AIC
estat ic
*878.2212

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln gdppc_ln gdppc_cng  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates
est sto gdppc4

*AIC
estat ic
*AIC=1329.814

****Create TABLE 20 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab gdppc1 gdppc2 gdppc3 gdppc4 using apx_gdppc.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 gdp_ln "GDP (ln)" gdppc_ln "GDP p/c (ln)" gdppc_cng "$\Delta$ GDP p/c" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Table 21: w/ homicide counts covariate

*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates
est sto count1

*AIC
estat ic
*AIC=1488.085

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
    y199* y2*   if durable2==0, cluster(ccode)

*store estimates
est sto count2

*AIC
estat ic
*AIC=129.943

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  hom_count  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto count3

*AIC
estat ic
*AIC=364.34
  
*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates
est sto count4

*AIC
estat ic
*AIC=1013.771

****Create TABLE 21 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab count1 count2 count3 count4 using apx_count.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 hom_count "Homicides (count)" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz
 
*Table 21: w/ homicide rates covariate
 
*Model 1: Full (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates
est sto rate1

*AIC
estat ic
*AIC=1527.282

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
    y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates
est sto rate2

*AIC
estat ic
*AIC=125.821 (but changes, for some reason)

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  hom_rate  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto rate3

*AIC
estat ic
*AIC=369.259

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates
est sto rate4
*AIC
estat ic
*AIC=1039.485

****Create TABLE 22 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab rate1 rate2 rate3 rate4 using apx_rate.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 hom_rate "Homicides (rate)" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

 ******************************
 *****IMPUTED MISSING DATA*****
 ******************************

*Get sense of how much missingness exists
 
sum seq_ln if polity!=. & inrange(year,1995,2014) & full==1
*n=3124; homicide data (n) drops it down to 1,982 (table 21) and 2,032 (table 22). 
* Over 1/3 lost.

*Impute Homicide Counts*

*Setup for imputation
mi set mlong

*Missingness table using model variables
mi misstable patterns confirmed seq_ln polity public_cor physical_vd ///
express_vd intensity2 info pop_ln hom_count  

*Tell Stata which variables you're imputing and which you're not imputing

*NOT IMPUTING
mi register regular confirmed seq_ln polity intensity2 durable2 pop_ln ///
 info express_vd physical_vd public_cor  y199* y2*

*IMPUTING
mi register imputed hom_count

*Imputation mode: use all right-hand side variables to impute  hom_count. 20 diff. datasets
mi impute mvn   hom_count ///
  = seq_ln polity intensity2  pop_ln info express_vd physical_vd public_cor, add(20) rseed(623)  force

*Summarize hom_count in different datasets (0= original, 1-20 are imputed datasets)
mi xeq 0 1 5 10 15 20: summarize hom_count   
 
*****************************
*Analysis using imputed data*
*****************************

*Set seed
set seed 45685219

*Table 23: w/ homicide counts covariate (imputed)

*Model 1: Full (NBREG)
mi estimate, post: nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates
est sto count_imp1

*Model 2: Auto (NBREG)
mi estimate, post:  nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
    y199* y2*   if durable2==0, cluster(ccode) nolog

*store estimates
est sto count_imp2

*Model 3: Ano (NBREG)
mi estimate, post: nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln  hom_count  ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto count_imp3

*Model 4: Demo (NBREG)
mi estimate, post:  nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_count  ///
  y199* y2* if durable2==2, cluster(ccode)  nolog

*store estimates
est sto count_imp4

****Create TABLE 23 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab count_imp1 count_imp2 count_imp3 count_imp4 using apx_count_imp.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 hom_count "Homicides (count)" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*Impute Homicide Rates*

*Need to clear dataset from previous imputation procedure
clear 
import delimited "Solis_fpa_july2019"
set more off

*Prepare dataset

*N. Log regime durability
gen seq_ln=ln(seq1)

*create dummy years
forvalues i=1992/2016{

gen y`i'=1 if year==`i'
replace y`i'=0 if y`i'==.
}
*

*Setup for imputation
mi set mlong

*Missingness table using model variables
mi misstable patterns confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  

*Tell Stata which variables you're imputing and which you're not imputing

*NOT IMPUTING
mi register regular confirmed seq_ln polity intensity2 durable2 pop_ln info express_vd physical_vd public_cor ///
  y199* y2*

*IMPUTING
mi register imputed   hom_rate

*Imputation mode: use all right-hand side variables to impute  hom_count. 20 diff. datasets
mi impute mvn   hom_rate ///
  = seq_ln polity intensity2  pop_ln info express_vd physical_vd public_cor, add(20) rseed(623)  force

*Summarize hom_count in different datasets (0= original, 1-20 are imputed datasets)
mi xeq 0 1 5 10 15 20: summarize hom_rate   
 
*****************************
*Analysis using imputed data*
*****************************

*Set seed
set seed 504504

*Table 24: w/ homicide rates covariate (imputed)

*Model 1: Full (NBREG)
mi estimate, post: nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
  y199* y2* , cluster(ccode)  nolog 

*store estimates
est sto rate_imp1

*Model 2: Auto (NBREG)
mi estimate, post:  nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
    y199* y2*   if durable2==0, cluster(ccode)

*store estimates
est sto rate_imp2

*Model 3: Ano (NBREG)
mi estimate, post: nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate ///
  y199* y2* if durable2==1, cluster(ccode)  nolog

*store estimates
est sto rate_imp3

*Model 4: Demo (NBREG)
mi estimate, post:  nbreg confirmed seq_ln polity public_cor physical_vd express_vd intensity2 info pop_ln hom_rate  ///
  y199* y2* if durable2==2 & inrange(year,1995,2014), cluster(ccode)  nolog

*store estimates
est sto rate_imp4

****Create TABLE 24 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab rate_imp1 rate_imp2 rate_imp3 rate_imp4 using apx_rate_imp.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" public_cor ///
"Public Sect. Cor., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)" fix "CPJ Unconfirmed" ///
 hom_rate "Homicides (rate)" ) ///
 varwidth(2) scalar(N_g) drop(y1* y2*) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

 *Table 25: w/ government compliance with judiciary covariate
 
*Need to clear dataset from previous imputation procedure
clear 
import delimited "Solis_fpa_july2019"
set more off

*N. Log regime durability
gen seq_ln=ln(seq1)

*create dummy years
forvalues i=1992/2016{

gen y`i'=1 if year==`i'
replace y`i'=0 if y`i'==.
}
*

*Model 1: Global (NBREG)
nbreg confirmed seq_ln polity comp_o physical_vd express_vd intensity2 info pop_ln ///
   y199* y2*, cluster(ccode) nolog 

*store estimates
est sto comp1

*AIC
estat ic
*AIC=2708.752 

*Model 2: Auto (NBREG)
nbreg confirmed seq_ln polity comp_o physical_vd express_vd intensity2 info pop_ln ///
   y199* y2* if durable2==0, cluster(ccode) nolog

*store estimates  
est sto comp2

*AIC
estat ic
*AIC=401.168

*Model 3: Ano (NBREG)
nbreg confirmed seq_ln polity comp_o physical_vd express_vd intensity2 info pop_ln ///
   y199* y2* if durable2==1, cluster(ccode) nolog

*store estimates  
est sto comp3

*AIC
estat ic
*AIC=979.4752

*Model 4: Demo (NBREG)
nbreg confirmed seq_ln polity comp_o physical_vd express_vd intensity2 info pop_ln ///
   y199* y2* if durable2==2, cluster(ccode) nolog

*store estimates  
est sto comp4

*AIC
estat ic
*AIC=1273.097

****Create TABLE 25 (Apx) for LaTex (basic table; I make some changes by hand once generated)****
esttab comp1 comp2 comp3 comp4 using jud_comp.tex, replace se aic obslast r2 ///
mtitle("Global" "Autocracy" "Anocracy" "Democracy"  ) ///
coeflabel(seq_ln "Regime-type Duration (ln)" polity2 "Polity Level" comp_o ///
"Court Compl., V-Dem" physical_vd "Physical Integrity, V-Dem"  ///
express_vd "Freedom of Exp., V-Dem" intensity2 "Armed Conflict" ///
 info "Information Flows" pop_ln "Population (ln)") ///
 varwidth(2) scalar(N_g) drop(y1* y2* _cons) b(%9.3f) t(%9.3f) r2(%9.2f) nolz

*****************************************
*****************************************
** SECTION G: COUNTRIES BY REGIME TYPE **
*****************************************
*****************************************

*Generate data for table 26

*Auto
tab country if durable2==0 & inrange(year,1992,2014)
*Ano
tab country if durable2==1 & inrange(year,1992,2014)
*Demo
tab country if durable2==2 & inrange(year,1992,2014)

*Finish log file
*translate Solis_fpa_Stata.smcl Solis_fpa_Stata.log, replace

****************
****************
**  END ~ jas **
****************
****************
